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DISCLAIMER 

The  findings  in  this  report  are  not  to  be  construed  as  an  official 
Department  of  the  Army  position  unless  so  designated  by  other  authorized 
documents . 

The  use  of  trade  name(s)  and/or  manufacturer(s)  does  not  constitute 
an  official  indorsement  or  approval. 


DESTRUCTION  NOTICE 

For  classified  documents,  follow  the  procedures  in  DoD  S200.22-M, 
Industrial  Security  Manual,  Section  11-19  or  DoD  5200. 1-R,  Information 
Security  Program  Regulation,  Chapter  IX. 

For  unclassified,  limited  documents,  destroy  by  any  method  that  will 
prevent  disclosure  of  contents  or  reconstruction  of  the  document. 

For  unclassified,  unlimited  documents,  destroy  when  the  report  is 
no  longer  needed.  Do  not  return  it  to  the  originator. 
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INTRODUCTION 
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i 

Shapiro  (ref  1)  states  in  his  introduction  that  he  does  not  intend  to  ! 

discuss  shape  preserving  properties  of  smoothing  or  iteration  of  the  smoothing 
process.  These  are  the  topics  which  will  be  discussed  here. 

The  basic  objective  in  this  report  is  to  take  a  piecewise  polynomial  of  low 
smoothness  (having  perhaps  no  derivatives)  and  replace  it  with  an  approximating 
function  having  any  desired  number  of  derivatives.  One  instance  in  which  one 
might  want  to  do  this  arises  from  computer-aided  manufacturing  where  one  might 
want  to  round  off  corners  in  piecewise  linear  geometries.  Another  instance 
arises  from  analysis  of  noisy  data  where  one  might  want  to  reliably  estimate  the 
second  derivative. 

We  will  adhere  to  the  idea  of  continuous  smoothing  through  integration  (as 
opposed  to  discrete  smoothing  through  summation)  because  it  becomes  a  trivial 
matter  to  interpolate  in  the  smoothed  function  or  any  derivative  thereof  even 
for  unequally  spaced  data. 

SHAPE  PRESERVATION  PROPERTIES 

Consider  the  averaging  operator  S  defined  by 

,  ,  1  ,x+h 

S  f(x)  «  --  /  f(t)dt  =  Fj (x) 

2h  x-h 

First,  the  operator  S  is  obviously  linear  because 

,  .  1  ,x+h 

S{af(x)  +  bg(x)(  =  --  /  af(t)  +  bg(t)dt 

2h  x-h 

1  ,x+h  1  ,x+h 

■  a  *  :r  /  .  f<t)dt  +  b  •  --  /  g(t)dt 

2h  x-h  2h  x-h 

=  aS { f ( x ) }  +  bS{g(x) } 

^Shapiro,  H.  S.,  Smoothing  and  Approximation  of  Functions,  Van  Nostrand  Reinhold 
Company,  New  York,  1969. 


Second,  S  preserves  1  and  x  because 


1  x+h  t  x+h  1 

S{1|  -  -r  /  ldt  =  --  |  -  --  (x+h-(x-h))  =  1 

2h  x-h  2h  x-h  2h 


and 


.  .  1  ,x+h  1  t!  .  x+h  1 

S{x}  «  --  /  .  tdt  -  --  --  |  «  --  ((x+h)2-(x-h)2) 

2h  x-h  2h  2  x-h  4h 


=  —  (x2+2hx+h2- (x2-2hx+h2 ) )  =  x 
4h 


Therefore,  S  preserves  all  linear  functions  because 

S{a+Bx}  =  S{A* 1+B*x}  =  AS{ 1 }  +  BS { x }  =  A+Bx 
This  is  the  extent  of  S's  accuracy  preserving  capabilities,  however, 
because  S  does  not  preserve  higher  powers  of  x  exactly: 


,  .  1  ,x+h  1  t3  x+h  1 

s{x2}  "  ::  /  u  t2dt  -  ::  I  •  ~  ((x+h)3-(x-h)3) 

2h  x-h  2h  3  x-h  6h 


1  h2 
—  (x3+3x2h+3xh2+h3-(x3-3x2h+3xh2-h3) )  =  xz  +  — 
6h  3 


Although  the  accuracy  preserving  ability  of  S  is  limited,  it  does  have  some  nice 
shape  preserving  properties  which  no  form  of  least  squares  approximation  has. 

For  instance,  S  preserves  monotonicity.  We  can  say,  for  example,  that: 

If  f  is  monotone  increasing  on  (A,B),  then  Fj  is  monotone  increasing  on 
(A+h,B-h).  Proof:  Assume  f  is  monotone  increasing  on  (A,B),  i.e., 


A<x<y<B  =>  f(x)  S  f(y) 
(implies) 


Let 


A  +  h<a$b<B-h 


By  definition, 


1  ,a+h  1  b+h 

Fl(a|  =  ii;  ;-h ,(t)dt "  is  Vh 


§ 
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U 
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In  the  second  integral. 


Therefore 


b-h  i  t  $  b+h 


A  <  a-h  $  b-h  $  t  <  b+h  <  B 


and 


We  therefore  have 


A  <  a-h  <  t  -  (b-a)  <  a+h  <  b+h  <  B 


A  <  t  -  (b-a)  <  t  <  B 

Since  f  is  monotone  increasing  on  (A,B),  we  also  have 

f (t-(b-a) )  <  f(t) 

The  second  integral  is  therefore  bounded  above  by 


b+h 

b-h 


f(t)dt 


but  this  is  just  F1(b).  We  finally  conclude  that 

Fl(a)  <  Fj(b) 

and  that  Fj  is  monotone  increasing  on  (A+h, B-h). 

S  also  preserves  convexity  or  concavity  "away  from"  inflection  points  in 
the  following  sense: 

If  f"(x)  exists  and  is  positive  for  A  <  x  <  B,  then  F^ (x)  >  0  for  A+h  <  x 

<  B-h.  Proof:  Assume  f"(x)  exists  and  f"(x)  >  0  for  A  <  x  <  B.  Take  arbitrary 

x  in  (A+h, B-h).  Therefore 

x-h  >  A  and  x+h  <  B 

Now  for  x-h<t<x+h,  we  have  A  <  t  <  B  and  f"(t)  >  0.  Therefore 

,x+h 

/  f"(t)dt  =  f ' (x+h)  -  f ’ (x-h)  >  0 


i  i1>  M  l^  U>'k*  «i*  ><*'#■* 


rf: 

s 


But  by  definition, 


,  .  1  ,x+h 

S{f(x)  =  F^x)  -  --  /  f(t)dt 
2h  x-h 


Using  Leibnitz's  rule  for  differentiating  integrals,  we  have 


Fi(x)  =  --  (f (x+h)-f (x-h) ) 
2h 


and  differentiating  again,  we  have 


Fl(x)  =  --  (f'(x+h)-f'(x-h)) 
2h 


F1(x)  >  0 


on  (A+h,B-h). 


S  also  preserves  positivity  in  the  following  sense: 

If  f(x)  >  0  on  (A,B) ,  then  S { f (x ) }  >  0  on  (A+h,B-h).  Proof:  Assume  f(x)  > 

0  on  (A,B).  Take  x  arbitrary  in  (A+h,B-h).  Therefore  x-h  >  A  and  x+h  <  B.  Now 
if  A  <  x-h  <  t  $  x+h  <  B,  we  have  f(t)  >  0  and  also  that 

,x+h 

/  f (t)dt  >  0 
x-h 


Therefore 


S { f  (x ) }  >  0 


on  (A+h,B-h).  A  corollary  to  this  theorem  is  that 

f (x)  >  g(x)  on  (A,B)  =>  S{f(x)}  >  s{g(x) } 
on  (A+h,B-h),  (S  is  a  monotone  operator.) 

Proof:  Assume  f(x)  >  g(x)  on  (A,B).  Therefore 

f(x)  -  g ( x )  >  0  on  (A,B) 


S{f(x)  -  g(x)}  >  0  on  (A+h,B-h) 
Linearity  of  S  gives  us  S { f (x ) }  >  S{g(x){. 


I'.ivi'iiV 


M.'ll.Wlt.'lU 


We  can  summarize  the  preceding  theorems  by  the  form: 

If  f(x)  has  property  P  on  (A,B),  then  S { f ( x ) }  has  property  P  on  (A+h,B-h). 
If  we  think  of  repeating  the  application  of  S  to  f  i  times  (S'* ) ,  we  can  easily 
prove  that 

f(x)  has  P  on  (A,B) 

*S> 

Si{f(x)}  has  P  on  (A+ih(B-ih) 

These  shape  preserving  properties  are  important  for  applications  in  industrial 
computer-aided  design  and  manufacturing,  and  in  data  analysis  situations  when 
one  or  more  derivatives  must  be  estimated. 


REPEATED  AVERAGING 

We  may  apply  the  smoothing  operator  S  repeatedly  in  order  to  obtain 
approximations  of  higher  smoothness  in  the  following  manner: 

,  .  1  ,x+h 

S{f(x)}  =  --  /  f(t)dt  *  F1(X) 

2h  x-h 


Fi+l (x)  »  S{F^ (x) }  i  >  1 

In  order  to  carry  out  this  process,  we  must  be  able  to  compute  the  successive 


indefinite  integrals  of  f 


f0(x)  a  f(x) 


f i+1 (x)  =  /  f i (t)dt  i  >  0 


We  now  compute  the  first  few  F's 


1  ,x+h  1 

Fl(X)  =  2h  ^x-h  f(t)dt  =  2h 


m 


u:u»A;iM4'a!i.n 


F2<x)  ■  i  l  C  Fi,,)dt  ■  £  C"  £ 


1  ,  x+h  ,x+h 

=  „  fi(*+h>dt  -  /  fi(t-h)dt 

(2h) *  x-h  x-h 


■  <f2(x+2h)  -  f2(x)  -  (f2(x)-f2(x-2h))} 

(2h) 2 


{f2(x+2h)  -  2f2(x)  +  f 2 ( x—  2h ) } 


f3(x)  * 


1  ,x+h  1  ,x+h  1 

—  /  F2(t)dt  =  —  /  -  ( f2 (t+2h)  -  2f 2 ( t )  +  f2{t-2h)  dt 

2h  x-h  z  2h  x-h  (2h) 2  *  z  z  f 


1  ,  ,x+h  ,x+h  ,x+h 

;:t;  U,-h  f2<'*2h>dt  - 2 *  /x.h  '2«-a»d‘f 


■-  1^3(x+3h)  -  f3(x+h)  -  2{f3(x+h)-f3(x-h) )  +  f3(x-h)  -  f3(x-3h)} 

I  3 


{f3(x+3h)  -  3f3(x+h)  +  3f3(x-h)  -  f3(x-3h) } 


The  appearance  of  the  binomial  coefficients  is  fairly  evident,  and  we  can 


guess  the  general  formula  for  the  ith  smooth  as 


Fi(x)  =  --™  J  (-l)k(^)fi(x+(i-2k)h/ 


All  we  need  to  do  to  prove  this  formula  in  general  (by  mathematical  induction) 


is  to  be  sure  that  it  is  true  for  i  =  1,  and  be  able  to  conclude  that  it  is  true 


for  i+l'on  the  assumption  that  it  is  true  for  i.  For  i  =1,  we  have 

1 


Fl(x)  =  J  (-l)k(J)fi(x+(l-2k)h)  =  ~  (fl(x+h)-fl(x-r 


which  we  have  already  shown  to  be  true. 


By  definition, 


1  ,x+h 

F<"(x>  ■  ih  U F«t)dt 


r 


wo 


mm 


Assuming  the  result  for  i. 


F<*i(x|  - ;;  Cl  1  (-Dk(’)f,(t.(-2kih)dt 


2h  x-h  (2h) 


(2h) i+l 


Tin  }  CH 


L„  \  (-l)k(^)  {f  i+1(x+h+(i-2k)h)  -  f  i+1(x-h+(i-2k)h)  j 


(2h)i+l 


Fi+l(x)  «  }  (-l)k(^){fi+l(x+(i+l-2k)h)  -  fi+1(x+(i-l-2k)h)} 


^  (-l)k(^)fi+1(x+(i+l-2k)h)  -  ^  (-l)k(^)fi+i(x+(i-l-2k)h)} 


(2h)i*l 


Ti;r  {  2  (_1)k<k)f i+l(x+( i+l~2k)h) 


(2h)i+l 


-  (-l)k"1(k]1)f i+i(x+(i-l-2(k-l) )h) } 


x  -h)i+l 


-t;-  {  }  (-l)k(k)fi+1(x+(i+l-2k)h)  +  }  (-l)k(k'1)fi+i(x+(i+l-2k)h)} 


At  this  point  we  must  mention  the  well-known  recursion  for  the  binomial 


coefficients 


<*-,>  *  <k>  ■  C> 


from  which  we  immediately  conclude  that 

<-!>  -  0  =  «ili> 

by  substituting  k  =  0  and  i+l  in  the  recursion. 


These  last  two  facts  enable  us  to  extend  the  ranges  of  summation  of  our  two 


! 


sums  appropriately 

i+1 

Fi+l(x)  =  (2h| i  +  1  {  2  (-l)k(^)i:i+l(x+(i+l-2k)h) 

k=0 

i+1 

+  ^  (-l)k(kl1)Fi+i(x+(i+l-2k)h)} 

k=0 

Therefore 

i+1 

Fi+l(x)  =  }  (-l)k{(k)  +  (kI1)}fi+l(x+(i+l-2k)h) 

k=0 

i  +  l 

k=0 


We  have  concluded  that  our  general  result  holds  for  i+1  if  it  holds  for  i,  and 
the  proof  is  complete. 

Having  established  the  validity  of  the  formula  for  the  ith  smooth  of  f, 

i 

Fi(x)  =  }  (-l)k(k)fi(x+(i-2k)h) 

k=0 

we  may  trivially  obtain  the  jth  derivative  of  Fn-  merely  by  subtracting  j  from 
the  subscript  of  f. 

i 

Fiix)  =  J  (-l)k(|J)fi_j(x+(i-2k)h) 

k=0 

If  f  is  not  differentiable  (and  usually  it  won’t  be),  we  must  be  sure  that  i  £ 

j- 


8 


a 

i 

Suppose  we  wanted  to  estimate  the  second  derivative  of  f.  We  would  of 

* 

course  need  j  =  2.  Suppose  in  addition,  we  wanted  the  second  derivative  to  be 

smooth  to  the  extent  of  being  differentiable.  The  lowest  value  of  i  that  we  \ 

could  use  would  depend  on  the  smoothness  of  f.  If  f  were  continuous,  then  fj  ! 

i 

would  be  differentiable  and  we  would  need  i  =  3  (at  least).  If  f  were  discon-  ' 

tinuous,  however,  f£  would  be  differentiable  and  we  would  need  i  =  4. 


SMOOTHING  ERROR  AND  CONVERGENCE 

In  this  section,  our  goal  is  to  obtain  a  general  result  describing  the  con¬ 
vergence  of  the  jth  derivative  of  the  ith  iterated  smooth  of  a  sufficiently 
smooth  function.  The  following  formula  for  integration-by-parts  will  be  used  a 
number  of  times: 

/  f(x)g(x)dx  =  f(b )/  g{x)dx  -  /b  f'(x)/X  g(t)dtdx 
a  a  a  a 

In  addition,  Leibnitz's  rule  for  differentiating  an  integral  containing  a  param¬ 
eter  will  subsequently  be  used. 


_d  ,£(x) 
dx  a(x) 


f (t,x)dt 


A 3(x)  3 

/  ,  ,  f(t,x)dt  +  f(0(x),x)0'{x) 
a(x)  dx 

-  f(a(x),x)a'(x) 


By  definition, 

S{f (x) }  =  -i  /X+h  f(t)dt 
x-h 


Letting  D  be  the  derivative  operator,  we  have  by  Leibnitz's  rule  for  differen¬ 
tiating  integrals: 


DSf(x)  =  (f (x+h)  -  f (x-h) ) 

but 

Of (x)  *  f'(x) 
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1  ,x+h 

SDf(x)  =  =£  /  ,  f’(t)dt 


=  2R  <*(x+h)  -  f  (x-h) ) 


*  DS  f{x) 

We  therefore  have  commutativity  of  the  S  and  D  operators  for  differentiable  f: 

SO  =  OS 

It  is  easy  to  show  that  in  general,  =  D-^S1  for  sufficiently  smooth  f. 

The  commutativity  of  the  smoothing  and  derivative  operators  tells  us  that 
when  we  apply  the  smoothing  operator  to  a  function,  we  simultaneously  apply  the 
same  operator  to  all  available  derivatives  of  the  function.  This  in  turn  tells 
us  that  the  shape  preserving  properties  of  the  smoothing  operator  extend  to  all 
derivatives  as  well. 


Recalling 


we  rewrite 


Jj)  1  ,x+h  (j) 

Fi+l(x)  *  —  /  Fi  (t)dt  j  >  o 
2h  x-h 


..(j)  ,x  (j)  ,x  (j) 

2hFi+1(x)  =  J  Fi  (t)dt  -  /  Fi  (t)dt 


Using  integration-by-parts  on  these  two  integrals,  we  have 

„  (j)  (j)  ,x  ,x  (j+1)  t 

2hFi+1(x)  =  Fi  (x)/  dt  -  /  Fi  (t)/  dudt 

x-h  x-h  x-h 


.  ( j )  ,x  x  (j+1)  t 

(Fi  (x)/  dt  -  /  Fi  ( t )/  dudt} 
x+h  x+h  x+h 


(j)  .x  (j+1) 

hFi  (x)  -  /  Fi  (t)  (t-x+h)dt 


(j)  ,x  (j+1) 

+  hFi  (x)  +  /  Fi  (t)  (t-x-h)dt 


Therefore,  we  get 


(j)  (j)  ,x  (j+1)  x  (j+1) 

2h(Fi+1(x)-Fi  (x))  »  -  /  Fi  (t)  (t-x+h)dt  +  /  F-j  (t)(t-x-h)dt 

x-h  x+h 

Using  integration-by-parts  again  on  each  integral,  we  have 

.  (j)  (j)  .  (j+1)  -x 

2h(Fi+1(x)-Fi  (x))  =  -{Fi  (x)/  t-x+h  dt 

x-h 

,x  (j+2)  ,t  .  (j+1)  ,x 

~  I  { t ) J  u-x+h  dudt)  +  Fi  (x)J  t-x-h  dt 

x-h  x-h  x+h 

,x  (j+2)  t 
-  /  Fi  (t )/  u-x-h  dudt 
x+h  x+h 

x 

„u/r.(j)  (j)  (j+1)  (t-x+h)* 

2h(Fi+1(x)-Fi  (x))  *  -Fi  (x)  - - -  w 

2  t=x-h 

t  x 

,x  (j+2)  (u-x+h)*  (j+1)  (t-x-h)* 

+  /  .  Fi  (t)  — - -  .  dt  +  Fi  (x) - - — 

x-h  2  u=x-h  2  t=x+h 

t 

.x  (j+2)  (u-x-h)* 

-  /  F*  (t)  - - - -  w  dt 

x+h  2  u=x+h 

h*  ^( j+1)  .x  (j+2)  (t-x+h)* 

*  ■  r F<  (x)  *  U F<  (t'  -  i  -  dt 

h*  (j+1)  .x  (j+2)  (t-x-h)* 

*  5-  F<  1x1  -  U  F’'  (t)  '  i  "  dt 

,x  (j+2)  (t-x+h)*  x+h  (j+2)  (t-x-h)* 

*  'x-h  F’  (t>  i  dt  *  K  Fi  (t)  “T“-  dt 

Using  the  absolute  value  triangle  inequalities  (for  sums  and  integrals),  we  have 

(j)  (j)  ,  ,x  (j+2)  (t-x+h)* 

2H  Fi+1(x)  -  Fi  (x)|  <  Jx_h  |  Fi  (t)| - - - dt 


,x+h  ,  (j+2)  .  (t-x-h)* 


We  must  now  define  the  following  function  norm: 

BgH(X;h)  =  least  upper  bound  of  the  set  {j  g(t)|  :x-h  <  t  <  x+hj 
=  max  {|  g(t)|  :x-h  <  t  <  x+h} 

if  g  is  continuous. 

Employing  this  norm,  the  previous  inequality  gives  us 

.  (j)  (j)  .  (j+2)  x  x+h 

2hj  Fi+1(x)  -  Fi  (x)|  <  %IIFi  #(x;h)(/  (t_x+h)*dt  +  /  (t-x-h)*dt) 

x-h  x 

(j)  h3  (-h)3  h3  (j+2) 

*  *,IFi  n(x;h)(- - - )  *  --  »Fi  « (X; h) 

Therefore,  we  have  another  preliminary  result 

.  (j)  (j)  .  h*  (j+2) 

lFi+l(x)-Fi  (x)h  --  "Fi  ll(X;h) 

The  special  case  of  i  ■  j  *  0  gives  us 

|  S{f  (x) }  -  f(x)|  <  y  Hf'  H  (x;h) 
o 

which  substantiates  our  previously  obtained  result  that  the  smoothing  approxima 
tion  is  exact  for  all  linear  functions.  Before  we  proceed  to  our  general 
results,  however,  we  need  just  two  more  preliminary  results.  The  first  norm 
theorem  is: 

If  ®fll(x;h)  <  HgH(x;k)'  then  HfH(X;h+l)  <  Hg»(x;k+1)  where  x  is  arbitrary 
and  h,  k,  1  >  0.  Proof:  Assume  flf n (X; h)  <  ltg#(X;k)-  By  definition, 

BfH(x;h+l)  s  max{j  f(t)|  :x-(h+l)  $  t  <  x+h+ij 
*  max{ Hf H (t;h) :x-l  <  t  <  x+l} 

<  max{llgB(t;kj  :x-i  <  t  <  x+ij 
=  max{|  g(t)|  :x-(k+l)  <  t  <  x+k+l} 
s  n9«(X;k+l) 


(note  that  use  of  the  triangle  inequality  does  not  make  this  inequality  par¬ 
ticularly  strong,  since  the  arguments  of  the  absolute  values  will  tend  to  have 
the  same  sign),  but  we  have  already  proved  that 


■  (j)  (j)  ,  h*  (j+2) 

I  p i +1  ( x )  “  Fi  (x)  I  *  “  ®Fi  **(x;h) 


Using  this  fact  in  the  previous  sum,  we  get 


|f5J;(x)  -  f(j)(x)|  <  ~  (IIF^i2)H(X;h)  +  IIF^2^»(X;h)  +  **-+^(j+2)ll(X;h)) 

Now,  applying  the  second  norm  theorem  to  the  first  term  i-1  times,  to  the  second 
term  i-2  times  etc.,  we  have 

I  F|J)  (x)  -  f(j)(x)|  <  ~  (Hf{J+2)B(x;ih)  +  nf(j+2)ll(x;(i-i)h)  +  ...-Hlf(j+2,"(x;h)) 

Since  the  first  term  is  obviously  the  largest,  we  have  our  final  result  for  suf¬ 
ficiently  smooth  f 

|fJj,(x)  -  f(j)(x)  I  «  ™  Iif(j+2)n(x;ih) 

O 

This  result  tells  us  that  the  derivatives  of  the  ith  smooth  converge  just  as 
fast  (quadratically  in  h)  as  the  ith  smooth  itself.  Also,  sines 

»fU+2)l'(X;ih)  <  llf(j+2)«(x;c0) 


(j+2) 


Uf ( j+2 ) j 


(X;oo) 


is  really  independent  of  x,  the  convergence  of  the  ith  smooth  and  all  its  deriv¬ 
atives  is  uniform.  Subsequently,  we  will  discover  that  this  error  bound  for  the 
jth  derivative  of  the  ith  smooth  is  the  best  possible  (smallest). 
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HIGHER  ACCURACY  SMOOTHING 

We  now  show  how  to  obtain  formulas  of  higher  accuracy  using  the  repeated 
smoothing  operator  S1 .  Unfortunately,  in  the  process  of  obtaining  higher 
accuracy,  we  must  give  up  the  shape  preserving  properties  of  our  approximations. 
Assume  in  what  follows  that  f(x)  is  a  polynomial  of  degree  n.  We  know  by 


definition  that 


,  ,  1  ,x+h 

S{f  (x) }  *  F^x)  =  --  /  f(t)dt 
2h  x-h 


Expanding  f ( t )  in  a  Taylor  series  around  x  gives  us 


1  ,x+h  r  f(k)(x) 


l  x+n  r  f'  'ixj 

Fi(x)  x  ~r  I  )  -------  (t-x)kdt 

2h  x-h  u  k! 


1  r  f(k)(x)  (t-x)k+1 


i  r  f  'ixj 

2h  k! 

k=0 


k+1  t=x-h 


1  r  f{K}  x  k+1  k+1 

=  --  )  - (hk+1-(-h)K+1) 

2h  ^  (k+1)! 
k=0 

1  v  f(k)(x)hk+1  k+1 

=  --  )  .  (1- (-1 )  1 ) 

2h  ^  (k+1)! 

k=0 

It  is  clear  that  every  other  element  of  this  sum  is  zero  (k  odd).  If  we  let 


G(x)  be  the  greatest  integer  <  x,  we  have 


Fi(x)  = 


k=0,2,4, . . . ,2G(n/2) 


f (k,(x)hk 

“"{k+l)T" 


G(n/2) 


f(2k)(x)h2k 
72k+l ) ! 


£ 

| 

Cjr 

is 


Letting 


and 


h2k 

{ik+i)T 


where  the  superscripts  on  the  C's  denote  the  smooth  index,  we  have 


Fl(x)  -  }  Ckf(2k)(x) 

k=0 


Since  F^  is  also  a  polynomial  of  degree  n,  as  well  as  all  successive  F's,  we  can 
write 


m 


k=*0 


and 


m 


fW(x>  *  }  d*1  fix?’  -  ] 


but 


k=0 


m-j 


j=0 


FS2j,,x, .  i  d  Air™  - 1  cj.jfi; 


k=0 


k=j 


i  + 


m 


k=0 


L 

j=n  ksi 


2  cjCk-jf(x 


Therefore 


•  I.I  I.*  I.' L«  I.1'!.*  I.'  li*  )!*,' 
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Carefully  interchanging  the  order  of  summation  in  this  double  sum,  we  have 


m  k 


Fi+l<x)  =  }  C^Mx)*  =  }  }  Cjcj-jfix)5 


k=0  j=0 


}  fix)*  ^  CjC^-j 


We  therefore  conclude  that 


*i+l  \  _i_i 
^k  *  2  ^j^k-j 


We  have 


co+1  3  CqCq 


ci+1  ■  Cqc{  +  cjcj 


c2+1  *  CqC2  ♦  CjCj  +  C2C0 


^m  3  C0Cm  +  CiCm_i  +  . . .  +  CmC0 


So  if  we  define  the  semi -convolution  product  of  vectors  a  and  b  as 


a*b  »  (a0laj,a2, . . .am)  •  (b0,bj,b2, . . .bm) 


(aobO'aobl  +  albO'aob2+albl+a2bO' • • 


aobm  +  albm-l  ♦  •••  +  ambo) 


we  can  write  C'i+1  *  C^C’,  where  removal  of  the  subscript  denotes  a  vector. 


We  use  the  term  "semi-convolution"  because  the  dimensionality  of  the  prod¬ 


uct  is  the  same  as  that  of  the  factors.  Also  note  that  Cg  _  j  for  an  ^ 


I 


I 


I 

» 


Let  us  now  derive  an  approximating  operator  which  is  exact  for  cubics 
f (x)  be  a  cubic  polynomial  (n=3).  Thus,  we  have  m  -  1.  We  can  write  the 
following  two  equations  and  subsequently  eliminate  f" (x): 


But  first. 


Fi  (x)  =  f  (x)  +  Cjf'fx) 
Fi+1(x)  =  f (x)  +  Ci+1f"(x) 

C1+1  *  (1,C}+1)  =  C’^C1' 

«  (l,Cj ) • (1 ,cj) 

=  ( i , cj+c| ) 


Therefore 


Cj+^  a  Cj  +  Cj 


ri+1  r1’  -  r1 
C1  *  C1  "  C1 

Summing  both  sides  of  this  equation,  we  have 


}  icj*1  -  cj>  •  }  c; 


cl+1  -  Cj  =  kcj 


Therefore 


We  have 


Cj+1  =  (k+1 )Cj 


_i  ,_i 
a  iCl 


Fi (X)  =  f (x)  +  iCif"(x) 


F1+1(X)  a  f (X)  ♦  (i+l)Cif"(X) 


-A  .w 


Eliminating  f"(x) 


Fi+l(x)  -  2?-  Fi(*>  ■  f(x)  -  f (x) 

iFi+1(x)  -  (i+l)Fi(x)  =  if (x>  -  (i+l)f (x) 

Therefore 


f(x)  =  (i+DF^x)  -  iF-j+j (x) 

=  (i+l)Sif (x)  -  iSi+1f(x) 

and  we  have  a  smoothing  operator  which  is  exact  for  all  cubic  polynomials 

s3  *  (i+l)Si  -  iS1+1 

We  will  now  obtain  an  operator  which  is  exact  for  all  quintic  polynomials. 

Since  n  =  5,  m  =  2,  and  we  write 

F-j  (x)  *  f  (x)  +  cjf”(x)  ♦ 

Fi+1{x)  =  f(x)  +  cj+1f*(x)  +  C3+1f |x| 

Fi+2(x)  *  f (x)  ♦  Cj+2f”(x)  +  C2+Zf [x| 

We  simply  eliminate  f"(x)  and  f[x|  from  these  equations,  leaving  f(x)  defined  in 
terms  of  the  three  smooths.  First,  however,  we  get  the  formulas  for  the  C's. 

Ci+1  -  (l,cJ+1,C2i+1)  =  C1 •C’ 

-  (l.CjjCo)* (i,c| ,C2) 

*  (i,cj  +  cj,c2  +  cjcj  +  C2) 

Therefore 

ci+1  *  c\  +  cj 

and 

-i +1  _ i  _ i _ i  _ i 

c2  *  C2  +  CjCj  +  c2 
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As  before, 


cl  *  ici 


_i+l 

C2  -  C2  =  i(Cj)2  +  C2 


Summing  both  sides,  we  have 


^  C2+1  -  C2  *  (c})*  )  i  +  )  Cj 


i=l  i=l 


Therefore 


„k+l  1  k(k+l)  <0i.,  .-1 

C2  -  C2  *  (Cj)*  +  kC2 


_k+l  k(k+l )  ,, 

C2  s  — - -  (Ci)*  +  (k+l)C2 


i  <M-1)  ,  ,  1 

c2  a  ------  (CJ)*  ♦  iC2 


Using  these  C's  to  eliminate  f" (x)  and  f 


ixj  ulti 


imately  gives  us 


f(x)  *  1j((i+l)(i+2)Fi(x)  -  2i(i+2)Fi+1(x)  +  i  (i+l)Fi+2(x) ) 

Our  smoothing  operator  which  preserves  quintics  is  therefore 

s5  *  3i((i+l)(i+2)S1  -  2i ( i+2)S1+1  ♦  i(i+l)S1+2) 

Using  s3  and  Sg  as  examples,  we  can  guess  that  the  operator 

k 

..i+k,  \  ,.k,  Si+J 

s2k+i  s  i(  k  )  )  (-i)J(i)  :■■■ 

*  *-  J  1+J 

j»0 

is  exact  for  all  polynomials  of  degree  2k+l  and  has  0(h2k+2)  error.  It  is 
interesting  to  note  that  these  s  operators  are  generalizations  of  Tukey's  (ref 
2)  method  of  twicing.  Tukey  uses  the  notation,  data  =  smooth  +  rough,  which  we 


2Tukey,  J.  W.,  Exploratory  Data  Analysis.  Add i son- Wes  ley.  New  York,  1977. 


»"A*j 


will  abbreviate  as  d  =  m  +  r.  He  then  defines  the  first  smooth  as  =  Sd,  and 
subsequent  iterated  smooths  by  m^+1  =  m-j  +  Sr.j,  where  r^  =  d  -  m-j .  Eliminating 
r.j,  we  get  m.j+i  =  Sd  +  m^  -  Sm^ .  For  twicing,  we  have 

m2  =  Sd  +  mj  -  Smj 
=  Sd  +  Sd  -  S2d 
=  { 2S-S2 )d 

For  thricing,  we  have 

m3  =  Sd  +  m2  ■  Sm2 

*  Sd  +  (2S-S2 )d  -  S(2S-S2)d 
=  (S+2S-S2-2S2+S3  )d 
=  (3S-3S*+S3 )d 

The  twicing  and  thricing  operators  are  the  same  as  S3  and  S5  with  i  =  1. 
Note,  however,  that  we  have  made  use  of  the  linearity  of  our  smoothing  operator 
S,  while  Tukey  (ref  2)  does  not  restrict  himself  to  linear  operators. 

KERNEL  FORM 

We  can  easily  write  the  first  smooth  in  kernel  form 


/  .  ,x+h  1 

S{f  (x) }  =  /  --  f(t)dt  =  /  Kj(t-x)f (t)dt 

x-h  2h  -oo 


where 


1 


-h  <  t  <  h 


2h 

0  elsewhere 


<i(t)  =  -j 

We  will  generalize  this  form  to  i  successive  smooths  and  find  K^(t)  where 

sMf(x)}  =  r  <i(t-x)f(t)dt 


We  will  need  the  following  simple  theorem  regarding  successive  integrations  of  f: 

^Tukey,  J.  W. ,  Exploratory  Oata  Analysis,  Addision-Wesley,  New  York,  1977. 
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Therefore 


s 1  f ( x )  =  )  C)  /x+(l'2k)h  f(t)  I*±li:?!<2!?ii22~  dt 

kto  <2h)"*  k  a  (1"1)- 

-  ^  /x+(l~2k)h  — i-l]k__  ((’)(x+(i_2k)h-t)i-lf(t)dt 

a  (2h) 1 (i-1) !  K 


,x+(i-2k)h 


gk(t)dt 


where 


K 

gk(t)  = - - - (')(x+(i-2k)h-t)i~lf(t) 

(2b)1 (i-1) !  K 

We  may  specify  a  arbitrarily,  as  long  as  it  is  independent  of  k.  Let  a  » 


ih,  hence 


r  x+(i-2k)h 

Snf(x)  *  /  /  9k(t)dt 

*•  v-ih 


l  r 

*-  x- 


x+(i-2k)h 


9k(t)dt 


S1f(x)  =  ^  /(l  2k)h  gk(t+x)dt 

u  -ih 


r  r ( i — 2k ) h 
2  /_.h  Gk(*)dt 


where  Gk ( t )  =  gi<(t+x). 


v.v.v-v  .v.v  - /.v'- 
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Writing  out  some  of  the  terms  of  this  sum,  we  have 


i  ,ih  ,(i-2)h  ,(i-4)h 

S'f (x)  =  /  G0(t)dt  +  J  ..  Gi(t)dt  +  /  G2(t)dt 

-in  -in  -in 


(i-6)h  -(i-2)h 

+  /  G3(t)dt  +...+/  G-j_i(t)dt 


-ih 


-ih 


Rewriting  this  sum  as  a  sum  of  integrals  over  disjoint  ranges,  we  have 

i  ,ih  ( i-2)h 

S  f(x)  =  ^«-  oik  Go(t)dt  +  /  G0(t)  +  Ga(t)dt 

( i-z } n  (i-4)h 


(i-4)h 

+  ( i-6)h  G°(t)  +  °l(t)  +  62(t,dt 


-(i-2)h 

+  •••  +  /  G0(t )  +  Gi(t)  +  ...  +  Gi-iftldt 

-ih 


Hence 


i-1  k 

i  T  ,(i-2k)h  r 

S1f (x)  *  )  /  )  G;(t)dt 

L  (i-2(k*l)  )h  ^  J 
k=0  j*0 


i-1  k 

X  ,(i-2k)h 

*  *  ( i-2(k+l) )h 
k*0  j=0 


]>  9j(t+x)dt 


=  2  /J’lrinh  2  — i"1-—  (j)  ((i-2j)h-t)i“lf(t+x)dt 
( 1-2 (k+1 ) )h  L  (2h) 1 (i-1) !  J 
k=0  j=0 

i-1  k 

}  f(t+x)  .  5  (-l)J’{t)((*l-2j)h-t),i“1dt 

4  ( i-2  (k+1)  )h  (2h)1(i-D!  L  J 

k=0  j«0 


but  we  want 


i  ,x+ih 

Snf(x)  =  /  K-j (t-x)f (t)dt 
x-ih 


ih 

=  /  ..  *i  (t)f  (t+x)dt 
-ih 
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We  can  therefore  conclude  that  for  0  <  k  <  i-1  and  (i-2(k+l))h  <  t  <  (i-2k)h, 


i(t)  =  -— * -  5  (-Dj(i)((i-2j)h-t)i“1 

(2h) 1 (i-1) !  *  J 


We  see  that  is  a  piecewise  polynomial  of  degree  i-1.  In  fact,  K-j  is  a  B- 
spline  (ref  3)  area  normalized  to  unity  with  constant  mesh  spacing  2h. 

Note  that  for  k  =  0  ( ( i -2 ) h  <  t  <  ih), 


Hence 


<i(t)  *  (ih-t) 


kSit)  oc  (ih-t)1"3"1 


( ih)  *  0  for  j  <  i-2 


ERROR  ANALYSIS  WITH  NOISE 


Let  DJ  denote  the  taking  of  the  jth  derivative  with  respect  to  x,  and  if 
denote  the  piecewise  linear  approximation  to  f  over  a  uniform  mesh  with  mesh 
width  t.  The  error  in  our  estimate  of  the  jth  derivative  of  f  is  given  by 

ej(x)  *  f(j)(x)  -  DVl(f+e)(x) 

where,  in  order  to  compute  our  estimate  of  f^,  we  take  the  sum  of  the 
underlying  function  f  and  noise  e,  evaluate  this  over  a  discrete  mesh  (take 
data),  define  a  piecewise  linear  approximation  to  the  data,  smooth  the  approxi¬ 
mation,  and  take  the  jth  derivative  of  the  smooth. 

We  may  rewrite  ej(x)  in  the  following  manner: 

ej(x)  *  f(j)(x)  -  DjSi(If(x)+ie(x)) 

»  f(j)(x)  -  D-Vlffx)  -  DVie(x) 

3de  Boor,  C. ,  A  Practical  Guide  to  Splines,  Springer-Verlag,  New  York,  1978. 
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but 

DjS1lf(x)  =  D-^S1  (If (x)-f (x)+f (x) )  *  DVf(x)  -  D-^S1  (f  (x)-if  (x)*) 
Therefore 

ej(x)  =  f(j)(x)  -  oVf(x)  +  O-Js^f (x)-lf(x))  -  D-jsiie(x) 

The  first  two  terms  denote  the  component  of  error  due  to  smoothing  alone.  Tne 
next  term  denotes  the  component  of  error  due  to  linear  interpolation  of  f,  and 
the  last  term  denotes  the  component  of  error  due  to  noise,  i.e.,  the  stochastic 
component.  We  denote  the  first  three  terms  by  A(x),  the  analytic  or  deter¬ 
ministic  component,  and  the  last  term  by  R(x),  the  random  or  stochastic  com¬ 
ponent  . 

Thus,  we  can  abbreviate 

ej (x)  =  A(x)  -  R(x) 

We  also  let  Ea,  Eq,  and  V  denote  the  arithmetic  mean,  quadratic  mean,  and 
variance  operators,  respectively. 

By  definition,  the  local  quadratic  mean  error  is  given  by 
Eq(ej(x))  =  (Ea(ej(x)*))*  *  qj(x) 

Therefore 

dj (x) *  *  Ea(A(x)*-2A(x)R (x)+R(x)«)  =  A(x)*  -  2A(x)Ea(R(x) )  +  Ea(R(x)*) 
but 

R(x)  =  D^s’ie(x) 
i  ,x+ 1 h 

*  DJ/  Ki (t-x)le(t)dt 
x-ih 

x+ih  ,  (j) 

=  /  (-l)JKi  (t-x)le(t)dt 
x-th 

(using  Leibnitz's  rule  and  the  fact  that  K^(±ih)  =  0  for  j  <  i-2). 
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Approximating  the  integral  by  a  sum  over  the  sampled  data  we  get: 


R(x)  ~  (-l)J  ^  KSj>(tk-x)e(tk)T 


The  expected  value  of  this  sum  is  clearly  zero,  and  using  elementary 
statistical  theory  for  the  variance  of  a  linear  combination  of  independent  ran¬ 
dom  variables  gives  us 

V(R(x) )  -  ^  (tk-x)2«r*T*  *  <7*t  ^  k\^  (tk-x)*r 


Re-approximating  the  sum  by  an  integral,  we  have 

,x+ih  (j)  ,ih  (j) 

V(R(x) )  -  a*T  j  K.-  (t-x)  *dt  «  <t*t  /  Ki  (t)  *dt 

x-ih  -ih 

This  approximation  will  be  good  for  r  small  relative  to  h. 

Since  Ea(R(x))  =  0,  Ea(R(x)*)  =  V(R(X)),  and  we  have 


where 


q j (x)*  *  A(x) *  +  V(R(x) ) 


.ih  (j) 

V(R(x) )  -  <t*t/  Ki  (t)»dt 

-ih 


and  o*  is  the  variance  of  the  noise. 


A(x)  =«  f(j)(x)  -  DJVf(x)  +  DJsi(f(x)-Jf(x)) 


and  taking  absolute  values  gives 


I  A(x)|  <  |  f(j)(x)  -  oVf(x)|  +  !  DV(f(x)-lf  (x)  )| 


The  previously  obtained  bound  on  the  first  term  is  given  by 


I  f<J,(x)  -  oVf(x)|<  ---  llf(j+2)ll(X;ih) 

6 


i 

$ 


We  now  obtain  a  bound  on  the  linear  interpolation  part 


|  D^S1  (f  (x)-lf  (x)  )|  *  |  0J’/  (t-x)  (f  (t)-lf  (t)  )dt| 

x-ih 


.x+ih  ( j)  x+ih.  ( j)  .  . 

|/  K*  (t-x)(f(t)-!f(t))dt|  <  /  .1^  (t-x)||  f(t)  -  if  {t)|  dt 

x-m  x-ih 


but  from  elementary  interpolation  theory,  we  know  that 


|  f  (t)  -  if(t)|  <  --  «f"B(t;T) 
o 


|  (f  (x)-!f  (x)  )|  *  ~  I  K^J 5  (t-x)!  Hf"8(t;T)dt 


*  ~~  ^  B (x; ih+r )  /  ..  I  Ki  (t-x)|  dt  =  --  lf"H  (X. -jh+r)  /\.  I  *i  (t)l  dt 


and  we  have  approximate  bounds  on  the  local  quadratic  mean  error  in  our  estimate 


of  the  jth  derivative  of  f: 


,iha  M+2i  T*  ,ih  .  (j) 

qj(x)a  <  {  afu  Z)n(X;ih)  +  --  if"ll(x.ih+T)J  |<i  (t)|dt}« 


-ih  (j) 

+  a*T^-ih  Ki  (t)*dt 


,ih  (j) 

I  ■  /  .u  I  Ki  (t)|  dt 

“In 


,ih  (j) 

J  -  /  K i  (t)*dt 
-ih 


and  approximate 


nf(j+2),l(x;ih)  ~  I  f(j+2)(x)|  +  ih  |  f(j+3)(x)| 
af,,|l(x;ih+T)  ~|f"(x)|+  (ih+T)|  f"'(x)| 


6S 


i 


u 


We  therefore  have  the  approximate  local  bound 


.ih*  ( j+2)  .  ( j+3) 

qj(x)*  <  {---  (|  f  (x)  |  +  ih  |  f  (x)  |  ) 

T*  I 

+  —  (I  f"(x)|  +  (ih+T)|  f"'(X)|  )  f  *  +  <7*TJ 
8 

and  the  approximate  global  bound 

qj a  r  qj(x)2dx  <  ■  /q  (■■■  (l  +  ih  I  ) 

J  L0J  L06 

t2I 

+  —  (j  f"(x)|  +  ( ih+T)|  f " r  (x)|  ) }  *dx  +  o*tJ 
8 

At  this  point,  we  will  consider  the  problem  of  estimating  I  and  J  (at  least 
for  j  *  0,1,2).  Obviously,  we  do  not  want  to  go  to  the  trouble  of  evaluating 
the  integrals  of  squares  or  absolute  values  of  derivatives  of  higher  order  B 
splines,  so  we  will  do  an  asymptotic  development.  First,  we  note  that  the  B 
spline  kernel  functions  (K^ )  appear  in  an  entirely  different  context  -  namely  in 
the  statistical  theory  of  obtaining  probability  distributions  of  sums  of  inde¬ 
pendent,  uniformly  distributed  random  variables  (ref  4).  If  xj,  X2,...x-j  have 
<1  as  their  common  probability  density  and  we  define  yi  by 


then  y.j  will  have  K-j  as  its  probability  density. 
We  can  easily  compute  the  variance  of  y^ 


V(Yi)  »  o-j*  =  ^  V(x|<)  =  i V ( x ) 


4Cramer,  H.,  Mathematical  Methods  of  Statistics.  Princeton  University  Press, 
New  Jersey,  1946. 
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but  x  has  density  Kj,  therefore  Ea(x)  =  0,  and 


We  then  have 


,h  x*  x*  .  h  h4 
V(x)  =  Ea(x4 )  =  /  —  dx  =  —  =  — 

a  -h  2h  6h  1  -h  3 


m  T 


We  now  bring  to  bear  the  central  limit  theorem  of  statistics,  which  tells 
us  that  the  probability  distribution  of  y,  will  approach  the  normal  distribution 
as  i  becomes  large.  In  practice,  this  normality  approximation  is  rather  good, 
even  for  small  i.  We  can  therefore  write 

K,(t)  - 

✓5? 

for  large  i.  First,  we  need  a  couple  of  derivatives  of  : 


Ki(t)  =  -  —  K-j  ( t) 


K'i'(t)  =  ---  ((-)4-l)K1(t) 
<7i4  oi 

The  incomplete  gamma  function  is  defined  by 


y(a,x)  *  /  ta_1e"tdt 
0 

We  will  subsequently  need  a  recursion  formula  for  y.  Using  integration-by¬ 
parts,  we  get 


$ 

& 
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P 

h 

\f 

t 


th* 

Pa  I 

ts* 
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H 
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y(a,x)  =  xa“l /  e-tdt  -  /  (a-l)t a-2/  e“ududt 
0  0  0 

a  xa-l  •  -e-t|  -  /  (a-l)ta_2  •  -e"u  |  dt 

0  0  0 

=  xa”l(l-e“x)  -  /  {a-l)ta-2(l-e-t)dt 
0 

*  xa-l(l-e_x)  -  /  (a-l)ta“2dt  +  /  (a-l)ta"2e-tdt 
0  0 

ta*l  j( 

=  xa~l ( i-e"x)  -  (a-1 )  -  |  +  (a-1)  /  ta_2e_tdt 

a-1  0  0 


=  xa"l(l-e"x)  -  xa“l  +  (a-l)y(a-l.x) 


Therefore 

y(a,x)  *  -xa~le"x  +  (a-l)y(a-l,x) 
We  also  have  the  ordinary  gamma  function 

T(a)  *  y(a,®) 

and  its  recursion 

T(a)  *  (a-l)r(a-l) 

Taking  j  *  0,  we  have 


I  -  /  I  <i(t)|  dt  -  /  Ki(t)dt  •  1  -  I0 
-00  -00 

J.f  K,(t).dt  .  /*  ----- 


— 00 

-oo  27T<7.j* 

,«  1 

,  1 

a® 

- e- 

■t*(7.dt  = - 

i  e“ 

-00 

2ttOi 

-00 

s  „L„  /27T 

.  ±  I 

-00 

27TOiy,2 

V2  it 

1 

iVn  o.j 
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K 


C5 


L*:« 
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Therefore 


where 


1  3*  Jo 

z/ii  i^h  i^h 


J0  =  ^3/n 


Taking  j  *  l,  we  have 


1  *  /  I  Ki(t)|  dt  =  r  ---  Ki(t)dt  *  r  ---  •  ------  e"^(t 

-  -"I*  -»  <V 

-<  —  *zr  -----  .-*faidt 

0  CTi  3V2it  0  oi*/2n 


cri/27t 


e^t*| 


o  1  , -  3*  _ _  Ii 

3  "  V2/n  *  —  •  Vl/H  =  — 
00  ffi  i^h  i%h 


(here 


Ij  =  /6/7T 


J  *  /  <i(t)*dt  =  C  ---  Ki(t)*dt  =  r  ---  •  -----  e'(t/< 

~°°  -^^i4  —  Oi*  2jrai* 


»  <7i*t*  1 


•  - - 


-00  (X.j 4  27T0.j* 


e”x  <7,-dt 


:T5  /_  t*e-t2dt 


1  1  i—  1  00 

I""  t«e-%t«dt 

2rrai 3  2V2  /27r  "® 


1  3 3 


4/ir  a.j*  a/H  i3/*h3  i3/^2f 


J1  3  t  /V* 


Finally,  taking  j  »  2,  we  have 


I  *  /“  I  K-(t)|  dt  =  / 


-oo  * 


,  t 

I  ( — )  * 

*i 


l|  K-j  (t)dt 


1 

<*ia 


t  . 
(")2-l| 
<?i 


-----  e~^t/a^2dt  =  if  ---  |  t *-l|  •  ------  e-^*oridt 

Vl jr  a,  0  ffi*  /irr  a. 


=  ---  Vlja  f  |  t*-l|  e-^t*clt  =  ---  JlJn  I/1  (l-t*)e-*t2dt  +  f  (t’-Ue-^dt) 

0^*0  0-j*  0  1 


Making  the  change  of  variable 

u  *  ijt* 
t  *  (2u)^ 

dt  ■  5i(2u)-,t2du  =  (2u)~^du 


I  *  ---  flTn(f  (l-2u)e-u(2u)_J<du  +  f  (2u-l)e-u(2u)-,idu} 

<Tia  0  H 

1  y,  n 

I  »  —  t/2/tt  {/  ( ( 2u ) -JS  -  (aujfye^du  +  /  ((2U)5*  -  (2u)-,*)e-udu} 

cr.j*  0  % 

a  ---  ✓27^{2-*t(*,*)  -  2^(V*,3i)  +  2*(r(V*)-y(s/*,l4))  -  2-*(rm  -  7(14,34))} 

Oi* 


a  —  y'2/7r{ -  2s/*y(V*,l4)  +  2*r(V*)  -  2~3ir(J4) } 

Oi* 

*  ---  /27^{2%y(34,l4)  -  23/i(-{34)*e-5*+>4y(54,34))  +  2*  •  34H14)  -  2“%r(34) } 

Oi* 

a  —  1/2 /n  •  2e_J4  a  —  fl/ne  =  2  •  —  ^2/Jte  =  — 
of-j*  o-j*  ih*  ih* 


where 


1 2  *  6/2/Tte 


J  *  /  K1"(t)*dt  *  /  —  ((--)*-l)*K1(t)*dt 


-oo  <7^* 


r  ---  ((--)2-l)2  •  -----  e‘(t/<Ti,4dt 

-oo  a.j  2itOi 2 

oo  i  i 

=  /  —  (t*-l)2  * - e^a-jdt 

-<*>  a  j*  2no.j2 


1  00 

-  /  (t*-l)2e-t*dt 

i*  0 


Making  the  change  of  variable 


t*  =  u 


t  = 


dt  a  Hu_5*du 

j  3  - - /  ( U— 1 ) 2 e~u  .  J$ij-^du  * - /  u~^(u2-2u+l)e“udu 


JNT-j*  o 


2;roi  *  0 


1  00 

- /  ( u3/2-2u^+u-5*  )e~udu 

2tto^  *  o 


Therefore 


1  ,  1 

J  3  - - -  (r(,/*)  -  2T( a/a )  +  r(J|))  -  .  (s/««*  T{\)-Z»\T(\)*T[X)) 

2  Jicr -j  *  2tto  -j  * 


r<35)  =  /n 


Therefore 


2na^* 


J  /  4  l/jr 


3  3*/» 


8/rr  04*  8/ri  i*/eh*  i*/*h* 


Summarizing  these  results  fcr  I  and  J,  we  have 


(h/T)J 


I0  =  1  ,  Ii  =  /iTw  ,  I2  =  6/2/ne 


(h/T)2J+1 


Jo  a  ,  Jj  =  »/aV3/tt  ,  J2  *  /37tt 
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Our  approximate  global  bound  on  qj  becomes 


1  ,L  ,ih*  (j+2).  .  ( j+3) 

Qj11  <  -  /  {-”  (I  f(x)  |  +  ih  |  f(x)  |  ) 
~  L  o  o 


t *  I j  a*rJj 

+ - (!  f”(x)|  ♦  (ih+r)|  f’”  (x)|  )}*dx  + 

8  ( h/i )  J  (h/iTJ  1 


Note  that  all  terms  in  the  preceding  bound  either  remain  finite  or  -  0  as  t 


-  0.  Also  note  that  for  a  given  t,  some  terms  ->  »  as  h  -  0.  What  we  would  now 


like  to  do  is  to  select  an  h  (in  terms  of  t)  which  roughly  minimizes  this  bound. 


We  would  then  like  all  terms  in  the  bound  to  -  0  as  t  (and  h)  -  0. 


Since  we  only  need  an  h  that  roughly  minimizes  the  bound,  we  can  freely 


neglect  small  terms.  We  can  begin  by  neglecting  the  f(j+3)  and  f”'  terms  as  h, 


t  -*  0.  We  then  have 


1  L  .ih*  (j+2)  T8Ij  .  .  .  «-‘-j 

qj* <  r  L  1  f(x)  1  +  -—/n 1  f,,(x>|  J2dx  + 

~  L  o  6  8 ( h/i ) J  (h/?)Zj+1 


ih*T*I 


T4I J  * 


-3"  j+2  +  "(h/T jj’  M2,j+2  +  64(h/7)23  42,2  +  (h/?)23;i 


where 


Hi,  j  *  [  /  I  fixlfliil  ^x 

L  0 


$ 

& 

& 


I  ^ 


35 


Therefore 


V  i ;;  ’ *h'^j+2, j+2  *  i\  <1'J/v"jT,Ij“2.j*2 

+  .1  T.h~2irj1  ,  +  h-(2j+l)i-(j^)a,Tj 

64  J  2,2  J 

=  ajh4  +  a2h2~^r*  +  a  h~2^T4  +  a  h~*2^+1*T 


where 


al  -  --  ^  *^j+2 , j+2 
1  l-i/2 

a2  ■  ;;  ’  Ij«2,j*2 


>3  <■  ~  <-JIj*»2,2 


a4  >  l-U+Wj.jj 


Differentiating  this  bound  with  respect  to  h  and  setting  the  result  equal  to 


zero  gives 

4a!h3  +  (2-j)a2h1-JT*  -  2ja3h_(2j+1)T4  -  (2j*l)a 4h_(2j+2)r  *  0 
Multiplying  through  by  h2J+2  gives 

4a1h2j+5  +  (2-j)a2hJ+3T*  -  2ja3hT4  -  (2j+l)a4T  =  0 
Clearly,  h  =  0  if  t  *  0,  but  we  would  like  an  h  in  terms  of  small  but  finite  t. 
It  seems  that  as  far  as  r  is  concerned,  we  should  be  able  to  neglect  the  a2  and 
a3  terms  of  the  last  equation.  We  will  now  do  this  and  subsequently  show  that 
it  is  justified.  Neglecting  the  second  and  third  terms  gives  us 

4a1h2J+5  -  ( 2j+l )a4T  =  0 

1 

51+5 


Therefore 


1 

(2j+l )a4  _ Jus 


luTVT 


■  ...  ...  *•«  «.«  ».•  ..»  «■»  ...  i.»  i..  u.mj 


where 


c ,  ,  ((2j*i)r»**1g‘j1)Si;S 

4al  1/9  i*Uj+2,j+2 

1  1 

=  _!  (9(2^1)^) 2j;g 

i  j+*/*fij+2<  j+2  /T  ^j+2,j+2 


Inserting  the  value  for  h  back  into  the  equation  in  which  we  neglected  the 
second  and  third  terms  gives  us 

j+3  1 

4a1c2J+5r  +  (2-j)a2T*cJ+3T  23+5  -  2ja3T4CT  2j+5  -  (2j+l)a4T  =  0 


5  j+13 


8j+21 


4aic2J+5T  +  (2-j)a2cJ+3r  23+5  -  2ja3cr  2J+5  -  (2j+l)a4T 


and  we  see  that  indeed,  the  a2  and  a3  terms  are  asymptotically  smaller  (as  t  - 
0)  than  the  two  terms  retained  (the  aj  and  a4  terms).  Also,  for  j  =  0,  the  a3 
term  is  zero  and  for  j  =  2,  the  a2  term  is  zero. 

We  now  insert  our  (roughly)  optimal  value  for  h  into  our  previous  bound 
expression.  The  bound  is 

qj*  <  ajh4  +  a2h2~JT*  +  a3h"2JT4  +  a4h‘(2J+1)r 

_4  2- j  -2j  -2 j— 1 

=  ajc4T  2j+**  +  a2T*c2_JT  2-i+®  +  a3T4c”2JT  2J+®  +  a4Tc"(2J+A)T  2-i+^ 


3  j+12 


6j+20 


3  ajc-x  2j+5  +  a2c2-JT2j+5  +  a3c_2jT23+5  +  a4C-(2j+D-r2j+5 

It  is  interesting  to  simplify  the  coefficients  of  the  powers  of  t  in  this  last 
bound  and  note  that  they  are  completely  independent  of  i  (the  amount  of 
smoothing) . 


9(2j+l )o* Jj  ---- 

k  =  ( . -)23+5 

^j+2, j+2 


£ 


v  v.v  >> 


mmm  w  m  j  w&w  j  p  i"  v  v  v  ^  ^  u.«  ^  v  r  w 
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Now 


ki — ^ . 


ale‘  '  j-  * *f*J+2,  j+2  •  k*<"' 
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a2C 


2‘J  =  -*  .  -I  uS-> 


24 


24 


J*  M2,j+2 


a3c"2j  *  --  i"jIj2ti2  2  *  *'2j^('2j)  s  '■  Ij*k‘2jy?  ? 

64  J  2'2  64  J  2'2 

a4c-(2J+1)  «  rlj^la.jj^^Dk-!2^)  « 

Our  roughly  optimal  asymptotic  error  bound  is 

4  3j+12 


*  <  ~  k4 


4j+2 


,j+2T2j+5  ♦  --  Ijk2~ju2, j 


6j+20 


j+2T 
4 


2j+5 


♦  --  Ii*k-2j42  2T2j+5  ♦  o*Jik-(2j+l)T2j+5 
64  J  '  J 

Clearly,  the  second  and  third  terms  play  a  minor  role  as  far  as  t  is 
We  can  therefore  further  simplify: 


qj*  <  (--  k4pj+2f  j+2  +  <7*Jjk-(2j+l))T2j+5 

**  v  O 


Since 


9(2j+l)a*Jj  ---- 

k  =  ( - -)2j+5 

yj+2,j+2 


TOvTTTv 


concerned. 


38 


we  have 


1  9(2j+l)<7*  Jj  »;§  9(2j+l)o2 Jj  -  ™ 

qj*  <  Uj+2.j+2(— . -)  J  +  a4Jj( . -”)  J  >*  J 

~  36  ^j+2,j+2  ^j+2,j+2 


2j+l  4  4  2 j+1  2 j+1 

*  "  W j+2 , j+2^+^ ( ^ ( 2 j+1  )<7*tJ j ) 33+s  +  (a*TJj)^+®nj+2,j+2^+®(9(2j+l))  2^+5 
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2  j+1 


Therefore 


«j'  <  ~  ("")4j*2,j.22J*5(9(2j*l)‘',TJj)2j*5 


2  j+1 


Recalling  the  expression  for  our  nearly  optical  h  value 

1  9(2j+l)o*Tjj  --l- 

h<  «  —  ( - -) Zj+J> 

/T  ^j+2, j+2 

and  defining  the  constants  2 

Mj  =  (9(2j+l ) Jj ) 2j+5 


and 


we  summarize 


"j  '  i  (“")<(S(2j*l)Jj)2j*5 


Mj  <7*r  51+5 

hj  =  --  ( . )ZJ  5 

/T  *ij+2,j+2 


2  j+1  2 

... 


qj  <  NjMj+2j+2  J  (»“t) 


where  M0  =  1.345,  M2  »  1.53,  M2  =  1.74,  and  N0  =  0.674,  N2  =  0.597,  and  N2 
0.679.  Note  that  -  0,  h^.  -  0,  and  hj/T  -  «  as  t  -  0. 
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We  now  apply  these  two  formulas  to  the  particular  case  of  a  sinusoidal 


function 


f(t)  =  A  sin  wt 


First,  we  can  easily  compute  that 

A*u2J+4 
yj+2,j+2  *  2 - 


We  then  have 


2j+4  2J*2  2 

q  4  Nj - j2j+5 

J  ~  J  2 


If  we  now  define  the  error-to-noise  ratio  e  to  be  qj/er  and  the  signal-to- 


noise  ratio  s  to  be  A/a,  we  get 


qj  <  N j ( - )4J+10(a*T)2J+5 


4j+2  4  2j+l  2 

*  N-a4J;*°  +  2^(---- _ j  4 j+l5T2j+5 

J  9 


Therefore 


,  2  j+4  li:L 

«  <  - ,4j*10T2j-S 


To  insure  that  e  is  bounded  above  by  some  specified  number  E,  it  is  sufficient 

that  2  j+1  2 

N  (?!".-!. )4jn6T2j;5  =  E 
J  2 


E  “5",  2 

T  *  (--)  2  r 


2  j+1 
■)“ 


h  .  "J  („?:i — ,53^5  .  "A  .  "A  j« 

/T  4j+2,j+2  Ki  A*u2  J+4  /i  s*urJ 


» 


r» 

s 


*  » 4^  r» 


RJJ 

KJ 


Inserting  the  previously  computed  value  of  t  into  this  expression  for  hj  gives 


2j+l  1 


nj  2  e  "5  2 

hj  =  ✓?  (s*w2>4(n])  (sI«2j+4)  ) 


4  \2j+5 


2  j+5  2 j+5  1 


--  (( - ~77)  5  (")  5  )5j+S  -  < - - - 


/?  s*w2J+4 


✓?  {S*«2j;i)  (Nj:> 


Therefore 


E  *  2  k 


isNj  o>2j+4 


For  a  sinusoidal  function,  we  may  therefore  estimate  the  needed  sampling 


interval  t  and  window  parameter  h,  given  the  desired  error-to-noise  ratio  E,  the 


signal-to-noise  ratio  s,  the  amount  of  smoothing  i,  the  frequency  w,  and  the 


derivative  index,  j. 


ERROR  ANALYSIS  WITHOUT  NOISE 


Recall  from  the  last  section  that  the  approximate  asymptotic  error  bound  on 


the  global  error  is  given  by 


qj*  <  ajh4  +  a2h2-^T*  +  a3h'2jT4  +  a4h"(2j+1)r 


and  that  setting  the  derivative  of  this  bound  with  respect  to  h  equal  to  zero 


yields  (implicitly)  the  optimal  window  parameter 


4a1h3  +  (2-j)a2h1_JT*  -  2ja3h- (Zj  +  1  )t‘  -  (2j+l)a4h'(2ji+2)T  =  0 


where 


1 

«  *  aUj+2 , j+2 


-  -- 


j**2,j+2 


Ri 

$•; 

w 


HMUXVWi 


I 


1 


.1  juTVXv^  <* /r  *  k  "> 


>-l  U  »«!  < 


1  -4 

a3  =  H  i  Jlj^2,2 

a4  *  i'(j+Ji)a*Jj 

In  the  presence  of  noise,  the  aj  and  a4  terms  turned  out  to  be  the  most 
significant,  but  without  noise,  the  a4  term  is  zero  and  we  are  left  with 

4ajh»  +  (2-j)a2h1-jT*  -  2ja3h_(2j+1)T4  *  0 


If  j  =  0,  we  have 


4ajh3  +  2a2hT*  *  0 


and  we  see  that  there  is  no  real  optimal  h. 


If  j  =  1,  we  have 


4a}h3  +  a2rl  -  2a3h“*T4  *  0 


4a1h*  +  a2T*h3  -  2a3T4  *  0 

There  are  no  dominant  terms  in  this  equation,  so  we  may  not  neglect  any. 
Solving  this  quadratic  for  h3  gives 


If  j  =  2,  we  have 


h3  »  —  {-a2  +  (a22+32aia3)5*} 
8a^ 


4a^h3  -  4a3h"*T4  »  0 


i^h*  =  a3T4 


Therefore 


a3  ’/• 

h  =  (--)  t* 
al 

Inserting  the  expressions  for  the  a's  gives 

9I2*H2,2  i/#  t 

h  *  ( . )  (r>* 

16u4,4  i 


and 


q2a  <  ..  (|l2f4+(|l2#2l»4,4)%)Tt 

We  see  that  even  without  noise,  there  are  optimal  window  parameters  and 
corresponding  error  bounds  for  derivatives. 


BEST  ANALYTIC  ERROR  BOUND 

In  this  section,  we  show  that  the  previously  obtained  upper  bound  on  the 
analytic  error  in  the  jth  derivative  of  the  ith  smooth  of  a  smooth  function  is 
the  best  possible,  i.e.,  the  smallest.  We  will  streamline  the  error  analysis  by 
making  use  of  the  kernel  functions  and  their  properties 

DJVf(x)  *  sVf(x)  -  sMxi  »  /X+lh  Ki(t-x)f|Jtidt 

x-ih 

,ih  (j)  0  (j)  o  (j) 

*  /  ..  <i(t)f(t+x)dt  =  /  <i(t)f  (t+x)dt  -  /.  K-j  (t)f  (t+x)dt 
-ih  -ih  ih 

Using  integration-by-parts,  we  have 

Djsiffx)  »  f(j)(x)/°  K-j { t ) dt  -  /°  f(j+1)(t+x)/t  K-j (u)dudt 

-ih  -ih  -ih 


Let 


-  {f(j)(x)/°  <i(t)dt  -  /°  f(J+1,(t+x)/t  K-j  (u)dudt } 

in  in  ih 


I (t)  *  /  Ki (u)du 
-ih 


Therefore 


D-JS^f  (x)  =  f(j)(x)/lh  Ki(t)dt  -  J°  f (j+1) (t+x)I(t)dt 

-ih  -ih 


♦  / 


0 

ih 


f(j+l) 


(t+x) (I(t)-I(ih) )dt 


D^S 1  f ( x )  -  f  jx  j  =  -  /°  f (^+x|l(t)dt  +  /°  fH+x)(I(t)-l)dt 

-ih  ih 

Again,  integrating-by-parts, 

DVffx)  -  fix)  a  -  {f«T1)/°  I (t)dt  -  /°  f j4+xj/t  I(u)dudt} 

-ih  -ih  -ih 

♦  fix)1’/0  (I(t)-l)dt  -  J°  f^+xj/*  (I(u)-l)dudt 

ih  ih  ih 

-  I(t)dt  +  /lh  ( 1-1  { t )  )dt }  +  /°.  f^txl/  t  I(u)dudt 


ih  j+2  ih 

+  /  f(t+x)  /  (1-I(u))dudt 
0  t 


I(t)  =  /  K.j  (u)du 
-ih 


therefore 


-t 

I("t)  *  /  Ki(u)du  *  -  /  K<(-u)du 
-ih  ih 

-ih 

1  J  *i(u)du  a  I ( i h )  -  I(t)  =  1  -  I (t) 


Therefore 


DJSif(x)  -  fUl  *  fU)1){-  /°  I(t)dt  +  /lH  I(-t)dt) 

-ih  0 

,0  (j+2)  t  ih  (j+2)  ih 

+  /  f(t+x)  /  I(u)dudt  +  /  f(t+x)  /  I(-u)dudt 
-ih  -ih  0  t 

fHFN  -/  I(t)dt  +  /  I(t)dt}  +  /  f[^+x|/  I(u)dudt 

-ih  -ih  -ih  -ih 

ih  (j+2)  -t 

+  /  f(t+x)/  I(u)dudt 


iv4**r»  I’* 


SiKySvS1 


/lh  t  I (u)dudt  *  -  /  lh  fi-t5x)/t  I (u)dudt 

0  -ih  0  -ih 

./°  fW.,  I{u)dudt 
-  i  h  -  i  h 

We  thus  have 

O-Vffx)  -  f(j)(x)  =  /°  f  (J+2)  (t+x)/t  I(u)dudt 

-ih  -ih 

+  /  f(J+2)(x-t)/  I(u)dudt 

-ih  -ih 

Taking  absolute  values  and  using  the  usual  norm,  we  get 

|D-iSif(x)  -  f(j)(x)|<  2!lf(j+2)||(x.ih)/°_h  /*  I(u)dudt 

=  2«f(j+2)fl(X;ih)/^_h  K-j  (s)dsdudt 

By  reversing  the  order  of  integration  in  the  triple  integral  a  couple  of  times 
we  get 

0  t  u  Os* 

*i  (s)dsdudt  =  /  =-  (s)ds 

-ih  -ih  -ih  -ih  2 

,ih  <*i*  ih* 

=  %  /  s*Ki(s)ds  =  ---  *  --- 
-ih  4  12 

We  finally  get 

|oVf(x)  -  f(J)(x)|<  2«f<J-2>*(Xi<h)  .  ~I  »  lfl-’*2)l(x;ih) 

and  we  see  that  no  better  bound  is  possible. 


■v.-vi"w>’v3«T> 


A  COMPUTATIONAL  CONSIDERATION 


To  repeat,  the  ith  smooth  of  function  f  is  given  by 


S1-f(x)  =  — ^  (-l)kC)fi(x+(i-2k)h) 

(2b)1  K 

k=0 


where 


(  u- 1 ) 1  A 

fi<u>  =  Ja  f<*>  (i-l)T"  dt 


Now,  for  a  given  value  of  x,  a  is  arbitrary,  but  if  u  is  considerably  different 


from  a,  the  integral  f-j(u)  can  be  quite  large.  This  can  lead  to  a  major  loss  of 


significance  through  round-off  error  in  the  smoothing  sum.  To  avoid  this,  we 


prevent  f-j(u)  from  becoming  too  large  by  selecting  a  equal  to  x. 


Hence,  we  compute 


f-i(u)  a  /  <j>i(t)dt 


where 


(u-t) 

♦i(t)  *  f(t)  ^nr 


Now,  assume  f(t)  to  be  defined  piecewise  on  some  x  mesh.  Let  $  x  i  xj+i  and 


xm  <  u  <  xm+1  where  1  i  m. 


(if  1  >  m,  fi(u)  »  -  0-j(t)dt) 


if  m  =  I,  f,-(u)  =  /  d>-;  ( t ) dt 

x 


*1+1  u 

if  m  a  1+1,  f.j(u)  a  /  <J>.j(t)dt  +  /  <#>T  ( t  )dt 

X  Xj+1 


if  m  >  1+2, 


,xl+l  xl+2 

*i(u)  3  /  <<>i  (t)dt  +  /  0-j(t)dt 

x  x1+1 


r  m  fu 

+  ...  +  /  <#>1  ( t  )dt  +  /  <fr-j(t)dt 

xm-l  xm 


,xi+l  r  ,xn+l  .u 

/  +  )  /  *i(t)dt  +  /  4>i(t)dt 


nei+1 


We  therefore  need  to  compute  integrals  of  the  form 


-p 

/  $-j(t)dt  where  xp  $  a,/S  <  xp+1 


In  the  special  case  when  f(t)  is  piecewise  linear  and  continuous,  we  have 
f(t)  =  yp  +  (yP+i  -  vP)  for  xp  <  t  <  Xp+1 


where 


f(t)  3  yp  +  qp(t-x-p) 


yp+i  "  Yp 

qp s  ^;r:_xp 


Therefore 


<Mt)  3  (Vp  +  qp(t-Xp)) 


,  .0  (u-t)1’1 

fa  *i(t)dt  =  (yp  +  qp(t-u+u-xp) )  <*t 


V.V  V^TlVA-.V/.  j-iiVM  J-x 
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